function [ RP MD RPo] = simulation( initialRPavg , initialRPstd , ns1 , nd1 , si , Ni , DISP)

%% The function is initialized and called from the script 'Main.m'
%% The function returns 3 variables:
% RP - The final risk perception of all agents in the chain
% MD - The signal of the message at each chain position
% RPo - The initial risk perception of all agents in the chain


% Initialisation of the agents' initial risk perception RPo, chosen from a
% random distribution with mean and standard deviation given in the
% function parameters.
RPo = initialRPavg + initialRPstd.*randn(1,Ni);

% Risk perception is defined in the interval [0 1]
RPo(RPo>1) = 1;
RPo(RPo<0) = 0;

% Initialisation of the structure where agents' risk perception is stored
RP = RPo;


% Initialisation of the probabilities that a positive or negative statement is added or removed at each transmission).
KD = RP ;
KS = (1 - RP);

% Initialisation of the number of positive and negative statements
NS(1) = ns1;
ND(1) = nd1;
N = ND+NS ;

% Initialisation of the signal of the message, stored in the variable MBIAS
MD(1) = ND/N;
MS(1) = NS/N;
MBIAS(1) = MD(1)-MS(1);


% Beginning of the simulation, with variable i describing the chain
% position

for i=1:(Ni-1)
    

   % STEP 1:  the agent at position i receives the message and updates his risk perception:

        % Normalization of the valence of the message, to make it in the interaval [0 1].
        normMB = (MBIAS(i)+1)/2 ;

        %The agent i updates his risk perception according to equation [1]
        %defined in the main text of the paper.
        RP(i) = RP(i) + si*(normMB - RP(i));

        %Update of the probabilities of message mutation
        KD(i) = RP(i);
        KS(i) = (1 - RP(i)) ;

    
    % STEP 2: the agent at position i biases the message according to KS and KD
        
        % Variables tempS and tempD temporarily store the current numbers of
        % positive and negative statements in the message
        tempS = NS(i);
        tempD = ND(i);
        
        
        % Randomly add new positive or negative statement to the message
        if rand<KD(i)
            tempD = tempD +1 ;
        end
        if rand<KS(i)
            tempS = tempS +1 ;
        end        
        
        % Randomly remove an existing positive or negative statement to the message
        if rand<KD(i)
            tempS = tempS -1 ;
            if tempS<0
                tempS=0;
            end
        end
        if rand<KS(i)
            tempD =tempD - 1 ;
            if tempD<0
                tempD=0;
            end
        end 
        
        % Update the number of positive and negative statements at the next
        % chain position, as well as the total number of statements.
        NS(i+1) = tempS;
        ND(i+1) = tempD;
        N = NS(i+1) + ND(i+1) ;

    % STEP 3: The agent i transmits the message to the next agent i+1
  
        MS(i+1) = NS(i+1)/N;
        MD(i+1) = ND(i+1)/N;
        MBIAS(i+1) = MD(i+1)-MS(i+1);
    
end

% Update simulation variables for the last chain position
i=Ni;
normMB = (MBIAS(i)+1)/2 ;
RP(i) = RP(i) + si*(-RP(i) +  normMB );
KD(i) = RP(i);
KS(i) = 1 - RP(i) ;


% End of the simulation

%% Display the final results in a new figure
if DISP==1

    figure
    suptitle('Simulation results')
    
    set(gcf, 'Position' , [ 624         514        1192         494]);
    subplot(1,2,1)
    plot(RPo , 'ko' )
    hold on
    plot(RP , 'ko-' , 'MarkerFaceColor' , 'g')
    plot(MD , 'ro-')
    title('Collective Dynamics')
    legend({'Initial RP of agents' , 'Final RP of agents' , 'Signal of the message' }, 'Location' , 'SouthWest')
    ylim([0 1])
    xlim([0 Ni])
    xlabel('Chain position');
    ylabel('Risk level');
  
    subplot(1,2,2)
    plot(ND , 'ro-')
    hold on
    plot(NS , 'bo-')
    title('Content of the message')
    ylim([0 25])
    xlim([0 Ni])
    xlabel('Chain position');
    ylabel('Number of statements');
    legend({'Number of negative statements' , 'Number of positive statements'}, 'Location' , 'SouthWest')
    
    
end



end
